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I. INTRODUCTION 

Sensors and measurement devices are affected by the 
presence or strength of physical effects that influence 
their dynamics in a detectable way. A proper statisti- 
cal treatment of measurement data is important when 
inferring results from complex experiments. With the 
growing use of quantum systems for high precision mea- 
surements, a whole research domain of quantum metrol- 
ogy has emerged. Limitations to measurement precision 
from quantum mechanical uncertainties have been inves- 
tigated and protocols to use measurements to optimally 
distinguish differently prepared quantum state have been 
developed, cf. [1-6]. 

The continuous observation of a quantum system in- 
volves leakage of information via coupling of the system 
to a suitable meter, and an archetypal example is that of 
measurements of photons emitted from a quantum light 
source. Laser spectroscopy thus involves the excitation 
of a quantum system, and detection of the fluorescence 
signal as function of laser frequency permits a fit, e.g., 
to a Lorcntzian distribution and thus yields information 
about the resonance frequency and linewidth. The res- 
onance curve, however, represents only a part of the ac- 
quired data, as it omits details concerning the temporal 
dynamics and noise properties of the detection signal. 
With the emergence of stochastic Schrddinger and mas- 
ter equations, which determine the quantum state con- 
ditioned on the full, noisy detection signals, it has been 
a natural next step to develop strategies to extract in- 
formation from continuously probed systems. Immedi- 
ate applications then concern sensing of the magnitude 
of perturbations acting on the system, such as the mag- 
netic field probed in atomic magnetometers [7, 8], and 
near field effects, e.g., from nuclear spins, probed by a 
single NV-center in diamond [9, 10]. Theoretical strate- 
gies have been proposed to continuously update param- 
eter estimates based on the acquired data in a Bayesian 
manner on equal footing with the conditioned quantum 
state of the probe system, [11, 12]. The latter method 
is particularly useful, if the system can be approximated 
by Gaussian states [13, 14]. 



The purpose of the present paper is to provide a for- 
mal link between some of the central ideas in classical 
estimation theory and stochastic master equations, and 
to identify efficient and systematic means to estimate un- 
known parameters from quantum measurement records. 
We will, in particular, discuss and demonstrate methods 
applicable in cases where the parameter space is too large 
to permit a recursive Bayesian update procedure. The 
methods are general, but for concreteness, we will con- 
sider light emitting quantum systems, and we will present 
explicit analyses and results for direct photon detection 
and for homodync detection of the emitted radiation. 

In photon counting, the measurement signal is a dis- 
crete process N t , characterized by click events at specific 
times, and the density matrix p t of the emitter condi- 
tioned on the detection signal until time t satisfies the 
non-linear filter equation 



dp t 



-i [H, p t ] - - {c f c, p t } + Tr(c t cp t )p t 



cp t C 



Tr(ctcp t ) 



Pt 



dt+ 
W, (1) 



where the differential measurement result dN t is a Pois- 
son increment, which is either (no click event) or 
1 (detector click event). For the special case of a 
two level atom with upper (lower) states \e(g)) and 
with upper state lifetime I/7, the conditional expec- 
tation E[dN t \N t ] = Tr(Jcp t )dt, where c = y/j\g)(e\. 
The time-evolution of p t is therefore piecewise continu- 
ous (when dN t = 0), but interrupted by jumps p t n- 
cptc^ I Tr {$ cp t ) at discrete times (when dN t = 1). 

It is also possible to perform field amplitude measure- 
ments by homodyne and heterodyne detection. The mea- 
surement signal Y t is then a continuous function of time, 
and, e.g., when homodyne detection is performed on the 
fluorescence emitted by the decaying two-level atom, the 
conditioned density matrix satisfies an ltd stochastic dif- 
ferential equation 



dp t = [-i [H, p] - {c f c, p) /2 + cpc f ] dt+ 

(M( Pt ) - Tr(M( Pt ))pt)(dY t - Tr(M(p t ))dt), 



(2) 
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where M. (p) = cp + pc 1 . The differential measurement 
result dY t satisfies dY t = Tr(M(p t ))dt + dW t , where dW t 
is a Wiener increment with zero mean and variance dt. 

In a generic experiment, all terms in the stochastic 
master equation can be parametrized by a vector of clas- 
sical variables 9 e R™, such as laser-atom dctunings, 
which may in turn depend on the unknown values of ex- 
ternally applied fields, decay rates, temperature, etc. In 
order to solve Eqs. (1, 2), candidate values for these pa- 
rameters need to be specified, and the goal of parameter 
estimation by continuous quantum measurements is to 
identify the best candidate values for the parameters 6, 
given the actual measurement record (N t or Y t ). 

In Sec. II, we review general parameter estimation 
concepts relevant to this manuscript: Bayesian inference, 
likelihood functions, and the Fisher information. In Sec. 
Ill, we show how likelihood functions and the Fisher in- 
formation can be efficiently obtained from the solution 
of the stochastic master equation of continuously mon- 
itored quantum systems. In Sec. IV, we introduce the 
Markov chain Monte Carlo method for efficient sampling 
of the likelihood function in large search spaces, and we 
give numerical examples which illustrate the application 
and the results of our methods. In Sec. V, we present a 
conclusion and outlook. 



II. BAYESIAN INFERENCE 

Our theory of estimation is based on Bayes rule, 
P(D\6)P{6) 



P(6\D) 



P(D) 



(3) 



where P(6\D) is the probability density of the parameters 
9, given the observed data D. Informally, this object con- 
tains all the information about the system parameters 9 
contained in the observed data D. From this distribution 
we can calculate any estimate of interest, including the 
mean value, the mode and quantiles. An important ad- 
vantage of calculating the full probability density P(9\D) 
is that it explicitly contains information about the uncer- 
tainty of the estimates. 

Bayes rule relates the conditional probability density 
to the probability of observing the data given the param- 
eters P{D\9) and the existing prior information about the 
parameters P{9). The difficulty in using Eq. (3) stems 
from the denominator being a weighted integral over 
all possible parameter values P(D) = J d9P{D\9)P{9). 
This integral is high-dimensional when several param- 
eters are estimated, and the integrand can vary many 
orders of magnitude. 

To determine P(9\D) in practice, we therefore need a 
method for calculating P(D\9) and a numerically efficient 
method of calculating P(D). 



A. Likelihood functions 

Since the data D is a measurement record, i.e., a func- 
tion of time, its probability density, or likelihood, P(D\9) 
is difficult to define. Apart from its use in the Bayesian 
update rule (3), it is common to maximize the likelihood 
with respect to the parameters 9 and thus to estimate 
the true value of the parameter by the value for which 
the likelihood for generating the data is highest. 

Instead of maximizing P(D\9) with respect to 9, one 
may maximize the value of any strictly increasing func- 
tion / of P(D\9). The logarithm is commonly used, and 
the resulting function is then denoted the log-likelihood 
function. 

It is also possible to divide P(D\9) with any strictly 
positive function P (D), without changing the location 
of the maximum with respect to 9. Thus any function 
f(P(D\9)/Po(D)) can be used to determine the maxi- 
mum likelihood estimate of 9. We will use the term like- 
lihood function for any function L(D\9) = P(D\9)/P (D) 
and logarithmic likelihood for l(D\9) = \ogL(D\9). The 
likelihood function L(D\9) can to a large extent be chosen 
to have a convenient form. 

Since P(D) = J d9P(D\9)P(9) 
P (D) J d6L{D\6)P{6) we can rewrite Eq. (3) as 



P(9\D) = 



P{D\9)P{9) 



P (D) J d9L{D\9)P{9) J d9L{D\9)P{9) ' 



(4) 



In the following we shall calculate the likelihood asso- 
ciated with continuously monitored light emitting quan- 
tum systems, where the function Pq{D) is the probability 
density for either a Poisson- or a Wiener-process. 



B. Fisher information 

A reasonable question to ask is, how accurate is the 
Bayesian estimate on average, and what is the funda- 
mental limit on how accurate it is possible to estimate 
91 

The answer to this question is given by the Fisher In- 
formation matrix. The Fisher Information matrix is de- 
fined in terms of the probability density for the data given 
some parameter 9 and P(D\9) as 



/ = E 



9 log P(D\6) Y 
89 J 



(5) 



where the expectation is over all possible realizations of 
the data D. The Cramer-Rao bound [15] states, that any 
estimator for 9, 9(D) has a variance larger than 1/I(9 ), 
where 9q is the true value of 9. 

If one uses a uniform prior, the Fisher Information of 
P(D\9) is the reciprocal of the width of P(9\D) averaged 
over the possible measurement records. If a non-uniform 
prior is included, the reciprocal width of P(9\D) is then, 
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qualitatively, the sum of the Fisher Information for 9 and 
the reciprocal width of the prior. 

As described in section II A, we can use any likelihood 
function in place of the conditional probability P(D\9). 
The same result holds for the Fisher information. That 
is, in Eq. (5) we can use any likelihood function L(D\6) 
instead of the probability density. 

For multiple variables, the Fisher information is 
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E 



d\og L(D\6) d\og L(D\9) 
89, dOj 

,dL(D\9) dL(D\6) 
89, 89, 



L(D[ 



(6) 
(7) 



where L(D\9) is a likelihood function for observing D 
given the parameters 9. The Cramer- Rao bound now 
states (for any unbiased estimator 9), that K[(6i — 9f)(9j — 
Of)} > (1(9)-%. 



III. CONTINUOUSLY MONITORED 
QUANTUM SYSTEMS 

The jump and diffusion quantum filter equations (1, 
2) are special cases of the general transformation of open 
quantum system density matrices subject to the random 
back action of measurements. If, at a given instant of 
time, a measurement occurs with outcome to € M, there 
is an effect-operator Q(m) associated with each outcome, 
so that the state, conditioned on the outcome to reads, 



p\m : 



(8) 



Tr(Qt( m )Q( m )p)' 
The probability (density) for observing the result to is 

p(m ) =Tr(fi t (m)fi(TO)p), (9) 
and the effect operators obey the relation 

J dmfl\m)fl(m) = 1. (10) 

The normalization factors in Eq. (8) are exactly the 
probabilities (9) to obtain the corresponding measure- 
ment outcome, and a similar probabilistic interpretation 
holds for the non-linear terms including the coefficients 
Tr(c^ cp t )pt and Tr(M(p t ))pt, in Eqs. (1, 2). This implies 
that if the stochastic density matrix equation is solved 
without incorporating the renormalization factors, the 
decreasing trace of p with time yields the likelihood for 
the actual detection record to occur. 

In the quantum jump master equation, the jump prob- 
ability, and hence the decrease in density matrix norm 
associated with a single jump is proportional to the du- 
ration of the infinitesimal time step dt chosen for the 
simulation. This causes an undesired and inconvenient 
dependence of the likelihood function on dt and on the 
number of click events, that we can, however, eliminate 
by a simple extension of the theory [2] similar to [16]. 



We introduce an arbitrary positive function po (to) and 
rescale the effect operators Q(m) — > Cl(m) / y/p (m) so 
that they now obey the modified normalization condition, 

/ dmpo{m)Stf(m)Sl(m)=l, (11) 
Jm 

Eq. (8) still holds, but the probability distribution for 
the different outcomes factors 



p(m) = Po(m) Tr(Q\m)£l(m) p) , 



(12) 



and we have the freedom to choose the un-normalized 
conditional states, 



p\m = fl(m)pfl'(m), 



(13) 



whose trace depends now explicitly on the chosen func- 
tion po(m). 

The expectation value, denoted by E, of any function 
/(to) is given by 

E[/(m)] = / dmp(m) f (m) 
Jm 

= [ dmp (m)Tr(p\m)f(m)=E [Tr(p\m)f(m)}, 
Jm 

where E is to be understood as the expectation with 
respect to the reference probability p a . In the follow- 
ing, we will suppress the dependence on the measure- 
ment outcomes and simply write p rather than p\m for 
the conditioned density matrix. 

The trace of the conditioned state is renormalized by 
a factor that depends on the specific detection record 
and which does not change its relative dependence on 
different values of the unknown parameters 9. It thus 
still serves as a likelihood function for the Bayesian de- 
termination of their values. Our scaling with the func- 
tion po(m) in Eq. (12) is indeed equivalent to the scal- 
ing allowed in the definition of the likelihood function, 
L(D\9) = P(D\9)/P Q (D), in Sec. II A. The "ostensible 
probability" p n [2] provides a reference measure po(m)dm 
on the set of measurement outcomes, and for our appli- 
cation it serves as a convenient unit for the effect op- 
erators fl(m). The relative entropy of p with respect 
to po is given directly in terms of Tr(p) as S(p\po) = 
E[log(Tr(p))]. 

Describing continuous measurements as the N — >• oo 
limit of a process of N times repeated measurements, it 
is natural to consider a general reference probability dis- 
tribution on M N , such that the probability of the mea- 
surement record factors, 

p(m x ,...m N ) =ptf f \m 1 ,...m N ) 

Tr(n(m N ) . . . Q,(mi)prt(mi) . . . itf{m N )), (14) 

where convenient reference probabilities p^\m) = 
p (mi) . . .po(m N ) for the jump-type and diffusion-like 
measurements will be given below. 
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A. Jump type equation 



B. Diffusion equation 



For the jump-type measurements there are for each 
small time interval dt two possible detector outcomes, 
dN t = and dN t = 1. We use our freedom to choose 
Po as the probability for a Poisson process with rate A, 
i.e. p n (dN t = 1) = Xdt and p {dN t = 0) = 1 - Xdt, 
and the correspondingly normalized measurement effect 
operators 



ft = 1 - iHdt - -(c f c - A)dt 



(15) 



The probability for a detector click is p (dN t = 
1) Tr(f2{f2i/9i) = Tr(c^ cp t )dt as expected from a rate pro- 
cess, and the expected number of events is K[dN t \N t ] = 
Ti(c^cp t )dt, while the reference expected value is 
E [dN t \N t ] = Xdt. 

The un-normalizcd conditional quantum state can be 
expressed as follows 



dpt 



-i [H,p t ] - - {Jc,p t } + Xp t 



dt 



+ dN t 



X 



Pt 



(16) 



while explicit normalization leads to Eq. (1). 

The dynamics of p t is governed by the Hamiltonian H 
and the operator c which in turn depend on the parame- 
ters 9. The likelihood of a specific sequence of detection 
events at times t\, . . . ijv < t is simply L(ti, . . . t;v|0) = 
Tr(pt), and Eq. (16) thus provides a differential equation 
for the likelihood function L t = Tr(p t ) 



dL t = (XL t - Tr(c'cp t ))dt + dN t 



Tr(cTcpt) 
A 



(17) 



where we have suppressed L t 's dependence on 9. 

The solutions p t of Eq. (1) and p t of Eq. (16) obey 
p t — Ti(p t )pt — L t pt which can be inserted in (17) to 
yield, 

dL t = (A — Tr(J cp t ))L t dt + dN t [A -1 Tr( C t C p t ) - l] L t . 

(18) 

This shows that even though the likelihood is formally 
defined by the trace of the un-normalized conditioned 
density matrix, L t can be calculated from the normalized 
state p t satisfying Eq. (1). 

For numerical purposes it is convenient to work with 
l t — log £4 which satisfies 



For diffusion type measurements, describing, e.g., ho- 
modyne detection of light, the set of outcomes in a small 
time interval dt is the real numbers. We will here use 
the probability of a Wiener increment dW t , i.e. a normal 
distribution with mean zero and variance dt as our refer- 
ence probability pff . The effect of observing a result dY t 
is 



1 



Sl(dY t ) = l- iHdt - -Scdt + cdY t 



(20) 



and the probability for observing a given value dY t is 

P (dY t ) = p^(dY t )(l + Tr((c + J)p t )dY t ). (21) 

We can calculate E[dY~ t |Y" t ] = E [(l + Tr((c + 
J)p t )dY t )dY t \Y t ] - Tr((c + J)p t )dt and E[dY t 2 \Y t ] = 
E [(l + Tr((c + J)pt)dYt)dY t 2 \Yt} = dt, which implies 



dY t = Tr((c + c i )p)dt + dW t , 



(22) 



where dW t is a Wiener increment with respect to the full 
probability distribution p (while dY t is a Wiener incre- 
ment with respect to po). 

The un-normalized stochastic differential equation be- 
comes 



dp t = 



1 



-i [H, p t ] - ^ { ctc > Pt} + cf>tc ] 



dt 

+ {cp t + pt^)dY u (23) 

and it leads to the likelihood L t — Tr(p t ) satisfying 

dL t = Tv(M(p t ))dY t . (24) 

As above, we can also express L t in terms of the nor- 
malized solution to Eq. (2), 

dL t = Tr(M(pt))L t dY u (25) 

and the log-likelihood It — log L t satisfies 

dl t - Tv(M(pt)){dY t - Ti(M(p t ))dt). (26) 

C. Fisher information 

Using Tr(p t ) as our likelihood function, we can apply 
Eq. (7) to calculate the Cramer-Rao bound for estimat- 
ing the unknown parameters in the system dynamics un- 
der both types of measurements. Define the matrices 



Tr(pt) 



(27) 



where the derivative is with respect to the i'th component 
of the vector of parameters 9. The expectation value of 
Tr(/9j) Tr(/9j) with respect to the probability distribution 
dl t = (A — Tv^ cp t ))dt + dN t log(Tr(c T cp t )/A). (19) p ( Le the actual pro bability for generating a trajectory) 
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Figure 1: (Color online) The upper panel shows the three components of the Bloch vector for a two-level atom 
subject to laser excitation with Rabi frequency 57 = 1.3 and detuning A = 1.43 (dimensionless units). The atomic 
inversion is represented by the lower (blue) curve. The atom decays with a rate 7 = 0.55, leading to the observation 
of quantum jumps at the instants indicated by vertical dashed lines and the transient Bloch vector dynamics. The 
lower panel shows the estimation of the detuning A, treated as an unknown variable with a probability distribution, 
which is updated in Bayesian manner, conditioned on the measurement record. See text. 



will then be the ij-component of the Fisher Information 
matrix for the continuously monitored quantum system. 
We can therefore evaluate the Fisher information matrix 
numerically by simulating the stochastic master equation 
a large number of times and determine the expectation 
value Eq. (7). 

In practice, for the jump-type measurement, this re- 
quires solution of Eq. (1) together with a simultane- 
ous evaluation of the matrices p\ , which can, in turn, be 
determined from the inhomogeneous jump type master 
equation 



dpi 



1 



[H,pl\--{Jc,pl}+Tr(Jcpi) P l 



dt 



-i[d l H, Pt ]~-{d^c),p\} 



dt 



+ dN t {cp\J + {d lC ) Pt J + cp t (d^) -pi), (28) 

where the stochastic term dN t takes the same value as in 
Eq. (1), and where the derivative of the Hamiltonian and 
damping terms with respect to 9i are assumed known. 
Similarly, for the diffusion-type measurement 

dpi = [-i [H, pi] - {ct c , p\) /2 + epjet] dt 

+ [-< [diH, pt] - {di(c1c),pt} /2 + (9 iC )p t ct + cp t (d^)] dt 

(M(pl) + (d l M)( Pt )-TT(M(p t ))pl)(dY t -TT(M(pt))dt), 

(29) 

where the Wiener increment dY t — Tr(Ai(p t ))dt = dWt 
takes the same value as in Eq. (2). 

The Fisher information provides an average quantifier 
of the asymptotic uncertainty in the estimation prob- 



lem. With Eqs. (1, 28) and Eqs. (2, 29) we have shown 
how the Fisher information can be calculated by simulat- 
ing many independent sequences of the stochastic mas- 
ter equation for the two different types of measurement. 
These simulations have to be carried out for the candi- 
date values of the parameters to yield the precision ex- 
pected for an estimate based on a typical experimental 
run. As illustrated by comparison of other such precision 
measures in [11], different measurement schemes have dif- 
ferent resolving power, and in future work, we plan to 
address these differences in more detail, e.g., by compar- 
ing the Fisher information derived for the jump-type and 
for different diffusion-type measurements. 

We also note, that if the field/meter degrees of free- 
dom could be left unmeasured, the full entangled density 
matrix of the quantum system and the quantized radi- 
ation field would depend on the unknown parameters. 
Thus the general quantum Cramer-Rao bound derived 
by Braunstein and Caves [1] to determine a parameter, 
encoded in a quantum state, yields the ultimate accu- 
racy with witch the parameters in our state dynamics 
can be inferred using any type of measurements. Iden- 
tifying that accuracy, and investigating how closely it 
is approached by quantum jump and quantum diffusion 
measurements of the emitted light presents an interesting 
challenge for further studies. 



IV. NUMERICAL INVESTIGATION 

In this section we will illustrate the theory outlined in 
the previous sections with a few characteristic examples. 
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One approach for investigating P(6\D) is to compute the 
likelihood function L(D\6) on a grid. Using such a calcu- 
lation, posterior expectation values of 9 can be calculated 
by numerical integration. A numerical maximization rou- 
tine can also be used to find the maximum of L(D\9) and 
thus provide a maximum likelihood estimate of the pa- 
rameters. The uncertainty is given by the Fisher informa- 
tion found by solution of the stochastic master equation 
with samples of simulated detection records. The pos- 
terior probability density may have many local maxima 
and it can be difficult to find the global maximum of 
L{D\9) using standard maximization techniques. 

If the parameter space is very large, more efficient 
methods for sampling the likelihood function exist. To 
sample a function with an un-normalized probability den- 
sity tt(x), one can apply a random process for the can- 
didate values in the form of a Markov chain, where 
the values jump in an appropriately chosen manner so 
that they attain the correct relative probabilities. The 
transition probability t(x\ — > x 2 ) must hence be cho- 
sen such that it asymptotically reproduces the relative 
probability density ir(x). The requirement for the tran- 
sition rule t is then that the only function that satisfies 
/ dxf(x)t(x — > x') = f(x') is proportional to our desired 
tt(x). A generic way to construct such a Markov chain is 
the Metropolis-Hastings algorithm [17, 18] which is used 
in many areas of science, and we provide a brief review 
of our application of the method. 

The basic idea is to compare the relative probability 
densities of a randomly chosen candidate value x 2 with 
the one of the current value x\. The value X2 may be 
chosen randomly or, more conveniently, according to a 
proposal chain q(x\ — > x 2 ), e.g., in the neighborhood 
of x\. A correct sampling of the probability density is 
obtained by accepting x 2 with the probability 

a(x 1 ,x 2 ) = mm 1, r , (30 

V 7r(xi)q(xi -> x 2 ) J 

and otherwise retaining the value X\. If the proposal 
chain is able to explore the entire parameter space this 
Markov chain will have n(x) as un-normalized stationary 
distribution. 

A nice feature of the Metropolis-Hastings sampling 
method, is that it uses only ratios between different ar- 
guments of the functions tt and q. This implies, that 
we can use the un-normalized probabilities n(x), and for 
our purpose, we can use the likelihood functions found 
by solving (18) or (25) with the parameter values 9 = x\ 
and 9 = x 2 ). 

In summary, to sample the posterior density for the 
estimated parameters P(9\D) for a continuous quantum 
measurement using Metropolis-Hastings we select a ran- 
dom 9 from the prior distribution P(9), and we proceed 
as follows: 

1. Determine candidate 9 C according to some proposal 
distribution q(9 — > 9 C ). 



2. Calculate the likelihood or, equivalently, the log- 
likelihood l T for the data until the final time T, 
using the candidate 9 C and Eqs. (18, 25) or Eqs. 
(19, 26) depending on the type of measurement. 

3. Calculate ot{9,9 c ) — min(l,exp(^ — It)q{9 c — > 
9)/q(9 — > 9 C ), where It is the log- likelihood for the 
previous parameter 9. 

4. Accept candidate with probability a(9,9 c ), other- 
wise keep 9. 

These steps are repeated a large number of times, and the 
parameters sampled are then representative and can be 
used for determination of any property of the distribution 
P(9\D). 

In the simulations presented below, the proposal dis- 
tribution q(9 — > 9 C ) was chosen as a multivariate nor- 
mal distribution centered at 9 with a variance selected to 
achieve a reasonable acceptance rate of 10% to 50%. 

Many techniques exist for investigating the conver- 
gence rate and the correlation length of the Markov chain 
generated by the above technique [18]. In the simple ex- 
amples studied in the present manuscript, the conver- 
gence rate and correlation length are readily identified, 
but a more careful analysis of these issues is necessary 
when applying the technique to an experimental situa- 
tion with many parameters and uncertainties. 

V. EXAMPLES 
A. Two-level atom 

Consider a coherently driven two-level atom that de- 
cays by spontaneous emission of photons. The atom is 
described by the Hamiltonian H = (fl/2)a x + (A/2)a z 
and by a jump operator c = ^o-~ , where 7 is the ef- 
fective decay rate, er = {a x 1 a y , a z ) is the vector of Pauli 
spin-matrices, and a~ denotes the Pauli lowering opera- 
tor. The measurement record is the times at which pho- 
tons are detected with a photodetector. 

The top part of Figure 1 shows an example trajec- 
tory, assuming known values 7 = 0.55, = 1.3 and 
A = 1.43 for the atomic and field parameters (in dimcn- 
sionless units, e.g., relative to the decay rate of another 
excited state in the same atom). The continuous curves 
show the components of the Bloch vector r = Tr(/?er), 
and they display continuous evolution disrupted at dis- 
crete times, where discontinuous quantum jumps of the 
state occur associated with the detector clicks. In the 
bottom part of Figure 1, we have assumed that 7 and 
f2 are known, and we evaluate the probability distribu- 
tion for the detuning parameter on a grid, assuming a 
prior normal distribution for A with a standard devia- 
tion a a = 1.0 and mean value ha = 2.0. 

The A-distribution is conditoned on the same detec- 
tion record as applied in the upper part of the Figure, and 
we observe how the no-click periods cause a continuous 
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Figure 2: (Color online) The left panels show 
histograms of Markov Chain Monte Carlo sampled 

distributions of the parameters, fi, A, 7 in our 
two-level atom model. The prior knowledge of the 
parameters assumes normal distributions, shown by the 
dashed lines, with mean values /iq = 2.0, /ia = 3.0, 
fi K = 1.0 and standard deviations ctq = 0.8, cta = 1-0 
and a K — 0.5. The right panels display the correlations 
between the different pairs of sampled parameters. 



change of the posterior density, while the discrete jumps 
are accompanied by more abrupt changes, until the dis- 
tribution is well converged. The importance of the use 
of the whole signal and not only the mean photodctcc- 
tion rate, is easily understood by the observation that 
following each quantum jump, the atomic density ma- 
trix describes a transient damped Rabi oscillation, and 
the temporal probability distribution for the subsequent 
jump event is periodically modulated. Since the period 
of the transient modulation depends explicitly on Q and 
A the actual occurrence of the next jump strongly favors 
(disfavors) certain values of A and causes the conditional 
increase (decrease) in the probability density at those val- 
ues. 



With a single unknown parameter, it is possible to 
compute the likelihood function on a fine grid, but if we 
pass to the larger parameter space of more unknown vari- 
ables, we have recourse to more advanced search meth- 
ods. In Fig. 2, we show the results of running the Markov 
chain Monte Carlo-algorithm on the trajectory in Fig. 1 
with all three parameters Q, A, 7 treated as unknown. 
We assume normal distributed priors, shown with the 
dashed lines in the left panels of Fig. 2, and the his- 
tograms show the values for the three parameters sam- 
pled by the Markov chain. Since the trajectory is quite 
short, there is not sufficient information to perfectly in- 
fer the values of the parameters, and the joint densities 
of pairs of variables in the right panels indicate that two 
islands of likely values of the set of parameters are not 
resolved by the measurements. 

We have compared the distribution of time differences 
between click event in the rather short detection record, 
shown in Fig. 1, with the expected transient Rabi oscilla- 
tion dynamics and we find that they are, indeed, compat- 
ible with the different values for the pair of parameters 
f2, A, occurring with comparable probabilities in the up- 
per right panel in Fig. 2. With a few more "lucky clicks" , 
however, the distribution will favor one choice, and due 
to the correlations between our estimates for all three 
parameters, they may then rather rapidly all converge to 
the correct values. 

We have also calculated the Fisher information matrix 
for a photon counting experiment by applying the simula- 
tion methods described above, and we obtain the results 
shown in Fig. 3. The Fisher information matrix was cal- 
culated by simulating the stochastic master equation and 
the associated equations for the p\ for different choices of 
the parameters (fi, A) e [-3/2,3/2] x [-3/2,3/2], while 
the decay rate was assumed to be known and equal to 
7 = 0.55 in our dimensionless units and the initial atomic 
state was unexcited. 

We recall, that the Fisher information, evaluated at the 
estimated values il, A gives the uncertainties of these two 
quantities as well as their covariance. It is not surpris- 
ing that the sensitivity of the photon detection method 
depends on the actual values of the parameters. The 
fact that fi and A enter the problem as coefficients on 
non-commuting spin components in the atomic Hamilto- 
nian suggests that spin uncertainty relations may result 
in limitations on their joint determination, see also [13]. 
Such a fundamental limitation may be reflected by the 
apparent anti-correlation of the occurrence of large and 
small values of the Fisher information matrix elements 
7n,fi and 7 A)A in Fig. 3. 

The relative entropy between the signal probability 
p and a Poisson reference distribution po, S(p\po) = 
E[logL t ], i.e., the p-expectation value of logL t , is shown 
in Fig. 3(d). The reference distribution p has been cho- 
sen as a Poisson process with a rate set by the stationary 
emission rate X st = H, 2 -y/(^ 2 + 4A 2 + 2S1 2 ). The rela- 
tive entropy is close to zero in large regions of the Q, A 
parameter space, indicating that the emission process is 
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Figure 3: (Color online) Fisher Information matrix 
components for photodetection of a decaying two-level 
atom with decay rate 7 = 0.55 initially in the unexcited 
state up to T — 40. The two upper panels show the 
diagonal elements Iq,q and I a, A and the lower left 
panel shows the off-diagonal element In, A of the Fisher 
information matrix. The lower right panel, shows the 

relative entropy between the signal probability 
distribution p and po, where po is a Poisson process of 
with the intensity of the stationary emission rate for the 
two-level atom A st = fl 2 j/(j 2 + 4A 2 + 2ft 2 ), and 
7 = 0.55, see text. 



not very different from a Poisson process. In the regions 
with |f2| > 0.5, A « 0, the dynamics deviate significantly 
from a Poisson process due to the Rabi oscillations in 
the photon waiting-time distribution, and the ensemble 
of trajectories have a higher entropy. 



B. Bi- modal two-level atom 

Imagine now a situation, where the two-level atom is 
not subject to dynamics with a fixed set of unknown pa- 
rameters, but it may jump randomly between two fixed 
sets of values. Such jumps may occur due to changes in 
a binary variable in the surrounding environment, e.g., 
the quantum states \a) and \b) of a nearby atom, spin or 
mesoscopic qubit degree of freedom, or due to the atom 
moving spatially between two different positions in a laser 
field configuration. We will assume these state changes 
are purely classical, i.e. we neglect all coherences be- 
tween the configurations or positions \a) and \b), and we 
assume that both the Rabi-frcqucncy, the detuning and 
the decay rate of the two-level atom have different values 
for the two states. 



Figure 4: (Color online) Simulated signal from a 
bimodal two-level atom, undergoing jumps and coherent 
evolution with two alternating sets of parameters, a and 
b. The solid black curve is the (binned) observed signal 
while the red dashed curve shows the mean expected 
counts for the atom subject to the current set of 
parameters. The values used for this trajectory are 
il a = 1.1, A a = 1.3, 7 a - 1.6, n b = 2.2, A b = 0.2, 
76 = 2.4 and transition rates W(a — > b) = 0.03 and 
W(b^a) = 0.08. 



We describe the system using a conditional master 
equation where we include the environmental states 
I a) and \b) of the atoms in a block-diagonal den- 
sity matrix, p — p a ® \a) (a\ + /9& ® \b) (b\, where 
Pa (pb) is the density matrix for the atom associ- 
ated with the environmental state a (b). The system 
Hamiltonian is H = ((Q a /2)a x + {A a /2)a z ) ® \a) (a\ + 
{(n b /2)<T x + (A b /2)cr 2 ) ® 1 6) (b\ and the effective photo- 
detection jump operator is 



(V7o» (a\ + VYb\b) (b\) 



The transitions between the two configurations are de- 
scribed by incoherent jumping rates W(a b) and 
W(b — > a) and corresponding jump operators J a _>b = 
^W{a^ 6)1® 1 6) (a\ and J b ^ a = y/W(b -> o)l® \a) (b\. 
The system is now equivalent to an enlarged quantum 
system, and it is fully described as a single quantum sys- 
tem by the formalism outlined above. 

We have used the parameters il a — 1.1, A Q = 1.3, j a — 
1.6, £lb — 2.2, Ab — 0.2, 7^ = 2.4 and (slow) transition 
rates W(a -4- b) = 0.03 and W(b -4- a) = 0.08 to simulate 
a typical detection record for the system. In Figure 4, the 
black solid line shows the time-binned observed signal 
for this record as a function of time. As the changes 
between the two sets of parameter occur at low rates, the 
photon counting permits efficient Bayesian determination 
of the classical states a and b along the same lines as 
[19]. The red curve shows the mean photon scattering 
rates evaluated for the current estimate of which set of 
parameters applies. 

Treating all rates and coupling strengths as unknown, 
the large number of unknown parameters makes a 
straightforward Bayesian estimation of their values very 
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Figure 5: (Color online) Marginal distributions for the 
eight unknown parameters in our bimodal two-level 

atomic system. All prior distributions were taken to be 

uniform on the shown intervals as indicated by the black 
dashed line. The estimation was based on the actual 

click events of the trajectory, partly shown in Figure 4. 



complicated. In Fig. 5 we show instead the outcome of 
the Markov Chain Monte Carlo sampling of the eight pos- 
sible parameters over the same measurement sequence as 
in Fig. 4. All values were assigned uniform prior proba- 
bility distributions on the intervals shown (dashed lines in 
the figures), and the histograms show the concentration 
of the values sampled on the actual, correct parameters. 



VI. DISCUSSION AND OUTLOOK 

In this paper we have presented a general method 
for inferring the values of parameters that govern the 
time-evolution of continuously monitored quantum sys- 
tems. The systems are described by stochastic master 
equations, and we have shown that the trace of the un- 
normalized density matrix can be interpreted and applied 
as a likelihood function in standard statistical methods 



for parameter inference. Explicit differential equations 
for the likelihood function in terms of the normalized den- 
sity matrix are exemplified for the case of photon count- 
ing (19) and homodyne photodetection (26). The differ- 
ential equations for the likelihood allows us to use nu- 
merically stable and efficient stochastic master equation 
simulations as input to a variety of standard statistical 
estimation algorithms, e.g. Markov Chain Monte Carlo 
and direct maximum likelihood estimation. Our identifi- 
cation of the conditioned density matrix dynamics with 
the likelihood function, in addition, leads to an efficient 
method (28), (29) to simulate the Fisher information as- 
sociated with any particular measurement scheme, and 
thus to evaluate the confidence of parameter estimation 
by continuous measurements. 

We presented our formalism for the case of photodetec- 
tion, and in our examples we assumed that all emitted ra- 
diation is detected. If there are unobserved decohcrence 
or loss processes and, e.g., loss of the radiation signal 
before the detection, averaging over these processes sim- 
ply contributes further (deterministic) dissipation terms 
of' the Lindblad form V[J](p) = - { JU,p) /2 + J pjt 
in the master equations (1,2). The likelihood equations 
(19,26), however, remain unchanged. 

A technical element in our formulation of the the- 
ory involves the introduction of a reference probability 
Po, imposing a degree of freedom in the normalization 
of the effect operators and the density matrix condi- 
tioned on the measurement signal. The introduction of 
the reference probability pg is mathematically equivalent 
to converting the set of measurement outcomes into a 
classical probability space with a reference probability 
measure Pq and the quantum measurement effect oper- 
ators induce a probability measure on M via the rela- 
tion P(A) = f A dP (m)Tr(p\m) for subsets Ac M. In 
mathematical terms the likelihood function Tr(p|r7i), dis- 
cussed in section II, can then be identified as the Radon- 
Nikodym derivative dP/dPo(m). This points to a fur- 
ther generalization by transforming the probability mea- 
sure P to some other measure P such that dP/ dP = Z t , 
where Z t is a martingale, and where trajectories gener- 
ated with the transformed probability measure P should 
be weighted by Z t to obtain ensemble averages. This may 
provide a useful technique to control the variance in nu- 
merically calculated ensemble averages and to simulate 
the master equation and the Fisher information matrix 
more efficiently. 

The relative entropy S(P\Pq) between the probability 
measures P and Pq is the P-expectation value of log L t , 
S(P\P ) = E[logi t ] = J dP\og(dP/dP ) and we note 
that the Fisher information is nothing but the relative 
entropy between Pg and P' g for infinitesimally close 9 and 
0' . Apart from their importance in parameter estimation, 
emphasized in this manuscript the entropy S(P\Po) and 
the Fisher information 7,j provide means to character- 
ize the stochastic dynamics of quantum trajectories in a 
manner similar to the use of entanglement susceptibility 
to characterize the correlations in quantum many-body 



physics [20-23]. 
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